library(tidyverse)
library(wesanderson)
library(MetBrewer)
library(sp)
library(sf)
library(classInt)
library(raster)
library(exactextractr)
library(fixest)

            
      #not all monitor data are publicly available so we do not release the raw monitor data only the processed correlations  

load("figdat/figsi9_data.RData")

pdf("figures/figure_SI9.pdf", width = 12, height = 12)

par(mfrow = c(3,3))
par(mar = c(4,5,4,3))
for(i in 1:length(wf2)){
  
  plot(data$wealth[data$country_w==wf2[i]], data$pm25[data$country_w==wf2[i]], pch = 21, col ='navy',cex=2, bg = data$bg1,  axes=F,xlab = "" ,ylab = "", ylim = c(min(c(data$pm25[data$country_w==wf2[i]], data$pm25_vd[data$country_w==wf2[i]]), na.rm = T)-1, max((c(data$pm25[data$country_w==wf2[i]], data$pm25_vd[data$country_w==wf2[i]])), na.rm = T)+1))
  
  
  points(data$wealth[data$country_w==wf2[i]], data$pm25_vd[data$country_w==wf2[i]], pch = 21, col ='red3',cex=2, bg = data$bg2,  axes=F,xlab = "" ,ylab = "")
  
  
  
  
  
  axis(1,cex=2)
  axis(2, las =2,cex=2)
  

  mtext(side = 1, text = "Wealth",cex=1, line=3)
  mtext(side = 2, text = "PM25",cex=1, line=3)
  abline( lm(pm25 ~ wealth,  data= data[data$country_w==wf2[i],], weight = 1/data$error[data$country_w==wf2[i]]), col = 'blue',lty=2)

  abline( lm(pm25_vd ~ wealth,  data= data[data$country_w==wf2[i],], weight = 1/data$error[data$country_w==wf2[i]]), col = 'red',lty=2)
  
    
  mtext(side = 3, text = unique(data$NAME[data$country_w==wf2[i]]),cex=2, adj =0)
 
 text(x = max(data$wealth[data$country_w==wf2[i]], na.rm = T)-0.1, y = max((c(data$pm25[data$country_w==wf2[i]], data$pm25_vd[data$country_w==wf2[i]])), na.rm = T), labels = paste("n=", sum(!is.na(data$pm25_vd[data$country_w==wf2[i]])), sep=""),cex=2 )  
  
}


dev.off()





